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Abstract 

We report numerical simulations of two-dimensional g-state Potts 
models with emphasis on a new quantity for the computation of spatial 
correlation lengths. This quantity is the cluster-diameter distribution 
function G^ ain (x), which measures the distribution of the diameter 
of stochastically defined cluster. Theoretically it is predicted to fall 
off exponentially for large diameter x, Gdi am <x exp(— x/£), where £ 
is the correlation length as usually defined through the large-distance 
behavior of two-point correlation functions. The results of our exten- 
sive Monte Carlo study in the disordered phase of the models with 
q = 10, 15, and 20 on large square lattices of size 300 x 300, 120 x 120, 
and 80 x 80, respectively, clearly confirm the theoretically predicted 
behavior. Moreover, using this observable we are able to verify an ex- 
act formula for the correlation length £d(fit) i n the disordered phase at 
the first-order transition point fit with an accuracy of about 1% — 2% 
for all considered values of q. This is a considerable improvement over 
estimates derived from the large-distance behavior of standard (pro- 
jected) two-point correlation functions, which are also discussed for 
comparison. 



PACS numbers: 05.50. +q, 75.10.Hk, 64.60. Cn, 11.15.Ha 



1 Introduction 



The physics of phase transitions is essentially governed by the behavior of the 
spatial correlation length £. While in some problems, e.g. at a continuous 
phase transition where £ diverges, it is often sufficient to know the qualita- 
tive behavior, there are also many applications which rely on quantitative 
estimates of £. This applies in particular to the finite-size scaling behavior 
near a first-order phase transition 1 where £ stays finite and sets the length 
scale above which asymptotic considerations should apply. 2 Since analytical 
predictions are scarce it is therefore of great practical importance to develop 
refined numerical methods for reliable computations of correlation lengths. 

In order to evaluate the accuracy of a newly proposed method one should 
apply it first to models where analytical predictions are available. The best 
known example is the two-dimensional (2D) Ising model where £ is exactly 
known at all temperatures. 3 But already the generalization to the 2D g-state 
Potts model 4 complicates the theoretical analysis considerably, and much 
less is known analytically. It was therefore a great success when a few years 
ago at least the correlation length £d(Pt) m the disordered phase at the first- 
order transition point f3 t for q > 5 could be calculated exactly. 5-7 Apart from 
heuristic arguments, no analytical predictions are available for the correlation 
length £ (Pt) in the ordered phase, and previous numerical simulations 8-10 
turned out to be difficult to interpret. This was the physical motivation to 
start a project 11 ' 12 with the goal to clarify conflicting conjectures for the ratio 
£,o/£,d at [3f The idea was, of course, to test the employed numerical methods 
first for the exactly known correlation length in the disordered phase 11 and 
then to proceed to the so far unexplored ordered phase. 12 

One often employed way to extract correlation lengths is to study the 
exponential decay of two-point correlation functions in the asymptotic limit 
of large distances. While this methods works perfectly for the 2D Ising 
model, for £d(Pt) °f the 2D g-state Potts models with q = 10, 15, and 20 
we experienced quite nasty systematic deviations from the exact answer by 
about 10% — 20%. 11 The deviations could be traced back to the unexpected 
importance of higher order excitations, but even though the Monte Carlo 
simulations were performed on quite large lattices and with a high statistics 
of about 50 000 — 100 000 uncorrelated measurements, least-square fits with 
sufficiently many correction terms turned out to be too unstable to predict 
reliable numbers. 
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A way out of this problem is to search for a different estimator of £ which 
is less affected by correction terms. A systematic search is certainly very dif- 
ficult, but one possible candidate was recently suggested in analytical work 13 
making extensively use of the Fortuin-Kasteleyn cluster representation 14 of 
the Potts model. In Ref. 13 it was shown that the distribution of the clus- 
ter diameter, Gdiam(^), decays exponentially for large diameter x, and that 
the decay constant is identical to the inverse correlation length (as defined 
from the decay of the two-point correlation function). This prompted us to 
investigate if the cluster-diameter distribution function is better suited for a 
numerical determination of the correlation length. In the following we report 
high-statistics Monte Carlo simulations of the models with q = 10, 15, and 
20, focussing on the properties of the new observable. As the main result 
it turns out to be indeed very well suited in the disordered phase, allowing 
for the first time a confirmation of the analytical formula for £d(Pt) with an 
accuracy of about 1% — 2%. Since we used larger lattices and considerably 
higher statistics than in our previous studies, 11 we discuss for comparison 
also the newly obtained estimates for £d(flt) from two different projections of 
the standard two-point correlation function. 

The remainder of the paper is organized as follows. In Sec. 2 we first 
recall the definition of the model and some exact results. We then discuss 
the simulation techniques and in particular describe the various estimators 
used to measure the correlation length. The results of our simulations are 
presented in Sec. 3, and in Sec. 4 we conclude with a brief summary of the 
main results and some final remarks. 

2 Model and observables 

In our Monte Carlo simulations we used the standard definition of the Potts 
model partition function, 4 

z = E e-/3£ ; E = - E <W s t = i,..., q , (i) 

M (ij) 

where j3 = J/ksT is the inverse temperature in natural units, i denote the 
lattice sites of a square lattice, (ij) are nearest-neighbor pairs, 5 StSj is the 
Kronecker delta symbol, and q is the number of states per spin. In all simu- 
lations we used periodic boundary conditions to minimize finite-size effects. 
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In the following we report results for the models with q = 10, 15, and 
20, employing lattices of size V = L x L with L = 300, 120, and 80, re- 
spectively. All simulations were performed in the canonical ensemble at the 
infinite volume first-order transition point f3 t = ln(l + y/q), at which the 
ordered and disordered phase can coexist. In a Monte Carlo simulation, the 
system can be biased into one of the two phases by the choice of the ini- 
tial spin configuration. To update the spins we used the Wolff single-cluster 
algorithm. 15 From a previous comparative study 11 we knew that in the disor- 
dered phase this algorithm clearly outperforms all other standard algorithms 
such as Metropolis, heat-bath, and Swendsen-Wang multiple cluster. 16 

The lattice sizes were chosen such that, for each value of q, L « 28^(A), 
with = 10.559519 . . ., 4.180954 . . ., and 2.695502 ... for g = 10, 15, and 

20, respectively. 5-7 Starting from a completely random configuration of spins 
it is then extremely probable that the system will stay in the disordered phase 
for a sufficiently long time, allowing statistically meaningful measurements of 
quantities being characteristic for the pure disordered phase. More precisely, 
by recalling that the escape probability V is proportional to exp(— 2a dL) 
and that in two dimensions the interface tension a ol i can be expressed in 
terms of the correlation length of the disordered phase, 7 a od = l/2£ rf (/? t ), 
one easily arrives at the order-of-magnitude estimate V oc exp(— L/^iPt)) ~ 
exp(— 28) ~ 10~ 12 . Finite-size corrections in the pure disordered phase are 
expected to be of the same order. 

In this work we mainly focussed on measurements of the probability dis- 
tribution of the cluster diameter, diamCj , which, in general, is defined as 
the maximal extension of a cluster in any of the D coordinate directions 
of a hypercubic lattice; for an illustration see Fig. 1. The cluster-diameter 
distribution function G iam (x) is then the probability, 

G dmm (a;) = /i(diamCj = x), (2) 

that the cluster Cj connected to a lattice site i$ has a given diameter x. 13 To 
increase the statistics we took advantage of the periodic boundary conditions 
and averaged G iam (:r) over all lattice sites i . In practice this amounts to 
recording a histogram H dl&m (x), whose entries at x — diamC are incremented 
by the size or weight |C| of each simulated cluster. If properly normalized, 
H dia,m (x) is then an estimator of the probability distribution G ia ™(x). As 
discussed in the Introduction the theoretically expected asymptotic behavior 



3 



of G diam (a;) in the disordered phase is an exponential decay governed by the 
correlation length ^, 13 

G diam (a;) = a exp(-x /&) + .... (3) 

By taking the logarithm of (7 diam an d performing linear two-parameter fits it 
is then straightforward to extract £d- 

For comparison we considered also in the new simulations the = 
2nn/L momentum projections (i = (i x ,i y )), 

g {n \i x ,j x ) = \Y, G (hjY^ ){iv - jy \ ( 4 ) 

with n = and 1 of the two-point correlation function 

G(i,j)^(5 SiSj -- q ). (5) 

For the measurements we actually decomposed the whole spin configuration 
into stochastic Swendsen-Wang (multiple) cluster and used the improved 
cluster estimator 17 

G(ij) = ^{e{ij)), (6) 

where 0(i, j) — 1, if % and j belong to the same cluster, and 6 = other- 
wise. In particular for small average single-cluster sizes (cp. Table 1), this 
procedure is more efficient than using directly the corresponding improved 
single-cluster estimator. 

As discussed previously, 11 to extract ^ from the large distance behavior 
of (4), non-linear four-parameter fits of the form 

gM(x) = g^(i x ,0) = ach(^^) + 6ch(c^^), (7) 

£d £d 

with 



e^er/Vl-(2vr<r/^) 2 - (8) 

are necessary. Below we shall report results for the first two projections with 
n = and n — 1. While the n = projection has been studied before on 
smaller lattices, 11 the use of the n = 1 projection in the disordered phase 
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is actually also new. Originally this projection was applied in the ordered 
phase where it is essential for removing constant background terms caused by 
the non-zero magnetization. 10 ' 12 Notice that for the large lattice sizes used 
in this study L m 28£ d , the difference in (8) between the fit parameter ffl 
and £d is only about 2.4%. 

The computer code was implemented on a T3D parallel computer in a 
trivial way by running 64 independent simulations in parallel. This allowed 
us to generate the very high statistics compiled in Table f . Here we followed 
the usual convention and defined V/(|C|)sc single-cluster steps as one Monte 
Carlo update sweep (MCS), where (|C|)sc is the average cluster size, and 
rescaled the integrated autocorrelation time of the internal energy, Ti ntiC , to 
this unit of time. Per r intiC we performed about two measurements of the 
projected correlation functions. The size and diameter of the clusters were 
measured for each generated cluster. The statistical error bars are estimated 
from the fluctuations among the 64 independent copies by using the standard 
jack-knife procedure. 18 The total running time of the simulations amounts 
to about five years of CPU time on a typical workstation. 

3 Results 

In all simulations we monitored the time evolution of the energy and magne- 
tization to convince ourselves that the system never escaped into the ordered 
phase. As a more quantitative measure we also computed energy and magne- 
tization moments which can be compared with exact 4 ' 19 or series expansion 
results. 20 The average and maximum cluster sizes and the maximum cluster 
diameter found in the simulations are given in Table 1. As a result of these 
tests we are convinced that, despite the very long run times, our results for 
£d can be identified with the pure disordered phase correlation length. 

The data for G iam (x) and g^°\x) are shown for q = 10, 15, and 20 
in the semi-log plots of Fig. 2. The continuous lines are one- and three- 
parameter fits to the Ansatz (3) and (7), respectively, with ^ held fixed at its 
theoretical value (= 10.559519 . . ., 4.180954 . . ., and 2.695502 ... for g = 10, 
15, and 20 5 ~ 7 ). Let us first concentrate on the new observable, the cluster- 
diameter probability distribution G iam (x). At first sight the constrained one- 
parameter fit to (7 diam looks less perfect than the constrained three-parameter 
fit to g(°\ since the data points are more randomly scattered around the fit. 
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The reason is that the correlations between the estimates at x and x + Ax 
are much smaller for G ^ n \x) than for g(°\x). This can be understood by 
noting that a cluster of diameter x contributes only to the one estimate of 
G iam (x) at x = x , but to all estimates of g^(x) with x < x (recall the 
cluster estimator (6)). 

The correlation length estimates resulting from various unconstrained 
two-parameter fits to G " am in intervals x min . . . x max with x max = 130, 50, and 
40 for q = 10, 15, and 20, respectively, are collected in Table 2. We see that 
the results are in very good agreement with the theoretically expected values, 
with only slight systematic deviations of about 1% — 2%. Contrary to the 
results 11 obtained from g^(x) the fitted values tend now to be overestimates 
for small x min . This tendency becomes obvious in Fig. 3 where we show the 
effective correlation lengths 

Q s (x) = l/HC(x)/C(x + l)}, (9) 

with C = (7 diam or g(°). The Q s (x) are just the inverse of the local slopes in 
Fig. 2. By recalling that neighboring values of G iam are much less correlated 
than those of g(°\ this explains the much larger error bars on the data for 
derived from G mm . Observe that Q s obtained from G mm develop a much 
more pronounced plateau for q — 15 and 20 than for q = 10, before also 
here the statistical errors increase and the data start to fluctuate around 
the theoretically expected value. To conclude this subsection, by using the 
cluster-diameter probability distribution as an estimator for the correlation 
length, we succeeded to confirm the theoretical prediction for £d(Pt) at a 
1% - 2% level. 

It would of course be unfair to compare the final estimates obtained 
from G iam (x) of the present study with the results from g(°\x) of previ- 
ous work 11 which used smaller lattices and lower statistics. In the present 
study we have therefore analyzed again g ( -°\x). Furthermore we discuss for 
the first time also g^'(x) in the disordered phase. In Fig. 2 we see that the 
constrained three-parameter fit to g^ yields an excellent description of the 
fall-off of g^ix) over more than four decades. Still, from an unconstrained 
four-parameter fit over the same x range with ^ as a free parameter we ob- 
tain for q = 10 an about 10% smaller value of ^ = 9.5(4). This confirms our 
earlier observation in Ref. 11 that four-parameter fits to g^ systematically 
underestimate ^. This is demonstrated in more detail in Table 2 where we 
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have collected the results of various fits in the intervals x m \ n . . . x max — L/2. 
For all three values of q we observe a clear tendency of increasing estimates 
for £ d with increasing x min . Still, the estimates in the last line for g^(x) are 
about 8% - 10% below the theoretical values. This tendency is also clearly 
visible in the behavior of the Q s (x) of shown in Fig. 3. Compared with 
Ref. 11 the statistics of the present simulations is higher by more than one 
order of magnitude. This allowed us to include larger x values in the fits 
and, as expected, improved the estimates of in particular for q = 15 and 
20. This clearly indicates that by further increasing the statistics also the 
remaining discrepancies could be removed. We can therefore conclude that 
there is nothing wrong, in principle, in using the standard two-point correla- 
tion function to estimate ^. Numerically, however, accurate estimates would 
require an enormous effort and would thus be a rather expensive enterprise. 

In Table 2 we also give the results of unconstrained four-parameter fits 
of the form (7) to g^\x) where, by using (8), we have already converted ^ 
to £d pa 1.02£^. We see that the estimates from g^ and g^> are strongly 
correlated, so that in the disordered phase nothing is gained by studying also 
the higher projections. 

Further investigations of the cluster-diameter distribution in the dis- 
ordered phase of the two-dimensional Ising and 3-state Potts models re- 
vealed, however, that the new observable is not always advantageous. Our 
results for the Ising model from a very long simulation of a 80 x 80 lat- 
tice with MCS/Ti nti e pa 12 288 000 = N meas in the disordered phase at f3 — 
0.70340888 pa 0.8(3 C are shown in Fig. 4. For (5 < f3 c the exact expression for 
the 2D Ising model correlation length is ^ = l/(/3* — /3), 3 where the dual 
inverse temperature j3* is given by (exp(/3) — l)(exp(/3* — 1) — q — 2. We 
see that here the ^ derived from G ' am clearly overshoot the exact value of 
£d = 2.6202906 . . . before they slowly approach it from above. Notice that 
P was adjusted such that agrees roughly with the value of the q = 20 
model at (3 t . The Q s of g^\ on the other hand, coincide with the exact value 
already for very small x, and a simple two-parameter fit of the form (7) with 
b = c = in the range x = 1, . . . , 40 = L/2 yields £ d = 2.62029(14), in 
perfect agreement with the exact result. A fit of G ' am according to (3) using 
only large x values in the interval x = 40, . . . , 56 gives a considerably higher 
estimate of ^ = 2.89(8) which, despite its large error bar, is only barely 
compatible with the theoretical value. The deviation is still about 10%. 
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If we choose (3 = 0.7891847 ~ 0.9(3 C (the dual inverse temperature of 
P* = 0.98), such that the exact correlation length is twice as large, ^ = 
5.2405812 . . ., we obtain qualitatively the same picture. This is shown in 
Fig. 5 for a simulation of a 160 x 160 lattice with MCS/r int)e « 3 200 000 = 
N me as- Here the linear fit of the data for g^ in the range x — 1, . . . , 80 = L/2 
yields ^ = 5.2400(6), again in nice agreement with the exact value. A fit 
of G iam in the interval x — 80, ... , 100 with £ d = 5.6(3), on the other hand, 
deviates again considerably by about 7%. 

Finally we show in Fig. 6 our results for the two-dimensional q = 3 Potts 
model at (5 = 0.95179503 « 0.95/? c (where 0* = 1.06). Here the lattice 
size is 160 x 160, MCS/r intje « 2 285 000, and A^ meas = 3 200 000. For small 
x we observe the influence of higher excitations in which, however, die 
out rapidly. Discarding therefore the smallest distances and choosing a fit 
interval of x — 7, . . . , 80 we obtain an estimate of ^ = 5.838(2), which is 
shown as the horizontal line in the lower plot of Fig. 6. Here the cluster- 
diameter distribution G iam is already slightly better behaved than for the 
two-dimensional Ising model, and a fit in the interval x — 80, ... , 106 gives 
a compatible value of ^ = 6.0(2), which now deviates only by 2% from the 
result of g(°\ 

4 Discussion 

Our numerical results clearly show that the cluster-diameter distribution 
G iam (x) is very well suited to extract the correlation length £d(Pt) of two- 
dimensional g-state Potts models with relatively large values of q. While 
analyses of the standard (projected) two-point function are plaqued by large 
systematic errors, with the new observable we succeeded for the first time to 
reproduce the theoretically expected values at a 1% — 2% level. 

For small values of q, however, the standard correlation function gives 
much more reliable results. For reasons not well understood to date, the two 
quite different correlators thus seem to behave complementary to each other. 

Also for the three-dimensional g-state Potts models with q — 3, 4, and 5, 
which undergo a first-order phase transition already for q > 3, our results 21 
for G iam and g^ in the disordered phase at the transition point (3 t as well as 
the corresponding effective correlation lengths look qualitatively as for the 
2D Ising model in Fig. 4. Also in these cases we found that g^ gives much 
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more reliable estimates of This suggests that the behavior of (7 diam does 
depend crucially on the value of q, but certainly not on the fact that the 
two-dimensional Potts models with q = 10, 15, and 20 were studied at their 
first-order transition point (3 t . The details of the 3D study will be published 
elsewhere. 21 
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Table 1: The average and maximum cluster size (\C\) and |C| max , the max- 
imum cluster diameter (diam C io ) max , the integrated autocorrelation time 
r int,c of the energy, the number of Monte Carlo update sweeps (MCS) in 
units of T inti e, and the number of measurements N meas . 



(3 — p t , disordered phase, 


single-cluster algorithm 




q = 10 


q = 15 


q = 20 




300 x 300 


120 x 120 


80 x 80 


(|C|) SC 


38.093(17) 


10.2254(11) 


5.87776(75) 


(|C|)sw 


2.726411(61) 


2.117866(18) 


1.854135(26) 


\C\ 

I ^ | max 


9353 


1908 


846 


(diam C k) ) max 


197 


73 


51 




w 59 


w 18 


~ 25 


MCS/r int , e 


600 000 


9 000 000 


4 200 000 


iV 

1 ' meas 


1600 000 


12 800 000 


10 395 000 
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Table 2: Results for £d(Pt) from two-parameter fits (3) to G diam (:r) for different 
fit intervals x min . . .x max , with x max = 130, 50, and 40 for q = 10, 15, and 
20, respectively. Also shown are results of four-parameter fits to g^(x) and 
g^(x) according to the Ansatz (7), (8), with x max = 150, 60, and 40 for 
q = 10, 15, and 20, respectively. 



observable 


q=!0 

^nn v 'inn 
oUU X ouu 


q = 15 
1 on v i on 


q = 20 

sn v sn 
oU X oU 


exact 


10.559519... 


4.180954... 


2.695502... 










^-^cliam 


40 10.90(2) 
48 10.92(3) 
56 10.88(3) 
64 10.84(5) 
72 10.75(8) 
80 10.6(2) 
88 10.5(2) 
96 10.3(3) 


20 4.297(4) 
23 4.286(5) 
26 4.267(8) 
29 4.26(2) 
32 4.25(2) 
35 4.27(3) 
38 4.25(4) 
40 4.23(6) 


13 2.766(3) 
15 2.761(3) 
17 2.752(5) 
19 2.744(7) 
21 2.74(1) 
23 2.73(2) 
25 2.70(3) 
27 2.68(4) 


gW(x) 


10 8.9(1) 
12 9.0(1) 
14 9.1(2) 
16 9.2(2) 
20 9.4(3) 
22 9.5(4) 
24 9.7(6) 


5 3.56(2) 

6 3.60(2) 

7 3.63(2) 

8 3.66(3) 

9 3.70(4) 

10 3.73(5) 

11 3.76(6) 


3 2.30(1) 

4 2.33(1) 

5 2.36(2) 

6 2.39(3) 

7 2.41(4) 

8 2.43(5) 

9 2.46(6) 


gW(x) 


10 8.88(7) 
12 8.96(9) 
14 9.1(1) 
16 9.2(2) 
20 9.3(3) 
22 9.5(3) 
24 9.6(4) 


5 3.55(1) 

6 3.59(2) 

7 3.62(2) 

8 3.65(3) 

9 3.69(3) 

10 3.72(4) 

11 3.75(5) 


3 2.293(7) 

4 2.33(1) 

5 2.36(2) 

6 2.38(2) 

7 2.40(3) 

8 2.42(4) 

9 2.44(6) 
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Figure 1: Illustration of the definition of the cluster diameter diamCj 
max{4,£ J/ }. 
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Figure 2: Semi-log plot of the cluster-diameter distribution G " am (:r) and the 
projected correlation function g(°\x) for q = 10, 15, and 20 at fit in the 
disordered phase. For clarity some data points are omitted for q — 10 and 
15. 
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Figure 3: Effective correlation lengths for q = 10, 15, and 20 at fit m 
disordered phase derived from the correlation functions shown in Fig. 2. 
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Figure 4: Correlation functions (upper plot) and effective correlation lengths 
(lower plot) for the 2D Ising model at /3 — 0.70340888. The horizontal line 
shows that exact value of &j = 2.6202906 
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Figure 5: Correlation functions (upper plot) and effective correlation lengths 
(lower plot) for the 2D Ising model at /3 — 0.78918147. The horizontal line 
shows that exact value of ^ = 5.2405812 
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Figure 6: Correlation functions (upper plot) and effective correlation lengths 
(lower plot) for the 2D 3-state Potts model at (5 — 0.95179503. The horizontal 
line shows our best numerical estimate of £ d = 5.838(2). 
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